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Ghost-free theories beyond the Horndeski class exhibit a partial breaking of the Vainshtein mech¬ 
anism inside non-relativistic sources of finite extent. We exploit this breaking to identify new and 
novel astrophysical probes of these theories. Non-relativistic objects feel a gravitational force that 
is weaker than that predicted by general relativity. The new equation of hydrostatic equilibrium 
is derived and solved to predict the modified behaviour of stars. It is found that main-sequence 
stars are dimmer and cooler than their general relativity counterparts but the red giant phase is 
largely indistinguishable. The rotation curves and lensing potential of Milky Way-like galaxies are 
calculated. The circular velocities are smaller than predicted by general relativity at fixed radius 
and the lensing mass is smaller than the dynamical mass. We discuss potential astrophysical probes 
of these theories and identify strong lensing as a particularly promising candidate. 


I. INTRODUCTION 

The acceleration of the cosmic expansion mm is one of 
the biggest mysteries in modern physics. For more than 
a decade, there has been a major theoretical and obser¬ 
vational effort (see [5] for a review) to understand the 
underlying mechanism. Infra-red modihcations of grav¬ 
ity are a promising and active area of research (see |4] for 
a review) but any viable modification must be compatible 
with solar system tests. This constraint usually renders 
any modification irrelevant on cosmological scales due 
to the need to tune any new effects to be sub-dominant 
to general relativity (GR) by five orders-of-magnitude or 
more. For this reason, theories that contain screening 
mechanisms (see [5] for a review) such as chameleon grav¬ 
ity 012], symmetrons [8] and the Vainshtein mechanism 
[5] are particularly promising due to their ability to give 
rise to order-one effects on cosmological scales whilst us¬ 
ing non-linear effects to recover GR in the solar system. 
The first two mechanisms (and their generalisations) pre¬ 
dict a plethora of astrophysical signatures that have been 
used to place incredibly strong constraints on the model 
parameters [1(1^20] . So strong in fact that modified grav¬ 
ity effects cannot drive the cosmic acceleration [2T] . This 
paper is concerned with the Vainshtein mechanism. 

Unlike chameleon-like theories, the Vainshtein mech¬ 
anism is so efficient at screening that, to date, there is 
an absence of novel small scale probe^ The Vainshtein 
mechanism features in many modihed theories of gravity, 
the quintessential paradigm being galileons [53]. Unlike 
chameleon-like theories, these are defined in the decou¬ 
pling limit and are the most general theory of a scalar 
field possessing the shift symmetry cj) ^ (p + -I- 6 in 

this limit. This symmetry leads to higher-order deriva¬ 
tive terms in the action such as (i9^()>)^n0 but ensures 
that the equations of motion are second-order, and hence 
the theory is ghost free. In order to examine the cosmol- 
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ogy of these theories, one requires a covariantisation so 
that the theory is dehned on any background. This co¬ 
variantisation is not unique and there are many inequiv¬ 
alent theories that give rise to the same theory in the 
decoupling limit. One covariantisation (if only the cubic 
galileon is present) is DGP gravity [53| where the galileon 
appears in the four-dimensional effective theory as the 
position of the 5D brane. This scenario is ruled out by 
cosmological observations [5^ and, furthermore, the self- 
accelerating branch has a ghost [5^. Another valid co¬ 
variantisation is ghost free massive gravity [27] and mas¬ 
sive bi-gravity |55| (see [l^ for a review of massive grav¬ 
ity). In its simplest form, massive gravity does not ad¬ 
mit exact Friedmann-Robertson-Walker (FRW) solutions 
[30] and more complicated forms do not have a standard 
cosmology [3T]. Massive bi-gravity does admit FRW so¬ 
lutions but these are unstable at the level of linear per¬ 
turbations |35H33] . One may instead try to covariantise 
the theory directly by promoting the partial derivatives 
to covariant ones. This process spoils the ghost-free na¬ 
ture of the theory since couplings to curvature tensors 
that are zero when the space-time is Minkowski give rise 
to higher-order derivatives of the scalar and the metric. 
One can remove these from the equations of motion by 
introducing non-trivial couplings of the scalar to curva¬ 
ture tensors in the action. One is then led to Covari¬ 
ant galileons [36] . Unfortunately, this covariantisation 
is ruled out by cosmological probes [57MD] . Covariant 
galileons are a subset of the Horndeski Lagrangian m, 
the most general Lagrangian with a scalar coupled to 
gravity that results in second-order equations of motion. 
Until recently, this was believed to be the most general 
ghost-free action for this particle content but [32 [33] (see 
also [33]) have shown that there are healthy extensions of 
this whereby the equations of motion are at first glance 
higher-order but reduce to second-order after the appli¬ 
cation of hidden constraints. In particular, the naive co¬ 
variantisation of galileons discussed above falls into this 
class of models. 

Interestingly, [45] have shown that the Vainshtein 
mechanism is only partially operational in theories be¬ 
yond Horndeski. Whereas it operates perfectly outside 
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sources there is an incomplete suppression of the fifth- 
force in the interior of extended sources due to the time- 
dependence of the cosmological field (see [3S1I17] for other 
examples of how cosmological time-dependence can be 
important on solar system scales in extended theories of 
gravity). This opens up the possibility of testing these 
theories on small scales using astrophysics. The aim of 
this paper is to discern the new novel probes that may 
potentially be used to probe these theories. 

We begin by looking at the effects on stellar structure. 
After finding the new equations for the metric poten¬ 
tials $ and inside sources, the equation of hydrostatic 
equilibrium predicted in these theories is derived. The 
theory introduces one additional dimensionless parame¬ 
ter, T, that depends on cj), the time-derivative of the cos¬ 
mological scalar and A, the mass-scale appearing in the 
quartic galileon Lagrangian. In order to gain some ana¬ 
lytic and physical insight into how the properties of stars 
are altered in these theories, polytropic models (with an 
equation of state P oc p^) are examined. These are useful 
because they decouple the effects of gravitational physics 
from those due to baryonic physics. The equations of 
stellar structure remain homology-invariant, which al¬ 
lows a Lane-Emden-like equation governing the structure 
of polytropic stars to be found. Main-sequence stars are 
well-described by 7 = 4/3 polytropes and so we pro¬ 
ceed to solve the Lane-Emden equation for this value. 
It is found that main-sequence stars are less luminous 
than their GR counterparts at fixed mass. In order to 
provide more accurate predictions that can be compared 
with data, we have modified the stellar evolution code 
MESA [33 HU to include the new equation of hydro¬ 
static equilibrium predicted by these theories. The code 
reveals that for low mass stars (~ IATq) the main effect 
of the modifications is that equivalent stars have higher 
effective temperatures at fixed luminosity. This effect is 
large on the main-sequence but is only minor on the red 
giant branch and is degenerate with keeping GR as the 
theory of gravity and increasing the metallicity. Higher 
mass stars exhibit similar features except the luminosity 
dimming is more pronounced. 

Finally, we turn our attention to galactic probes. We 
compute the modified rotation curves and lensing poten¬ 
tials for Milky Way-like galaxies. It is found that non- 
relativistic objects have smaller circular velocities at fixed 
radii than predicted by GR and that light is lensed less 
strongly. We comment on the use of strong lensing as a 
potential probe of this effect. 

This paper is organized as follows: Quartic galileons 
and the Vainshtein mechanism are reviewed in the next 
section where the equations governing the metric poten¬ 
tials $ and 4* are derived. In section |TTT| the new equation 
of hydrostatic equilibrium is derived. It is used to find 
the modified Lane-Emden equation governing the struc¬ 
ture of polytropic stars and we subsequently specialise to 
stars with equations of state of the form P oc i.e. 
main-sequence stars. The resultant equations are solved 
to show that main-sequence polytropic stars are dimmer 


than their GR counterparts. In the remainder of the sec¬ 
tion, the stellar structure code MESA is used to produce 
realistic predictions and it is found that in addition to be¬ 
ing dimmer, main-sequence stars are cooler (in the sense 
that they have smaller effective temperatures). Galac¬ 
tic effects are studied in section jT^ where we examine 
the rotation curves and lensing potentials of Milky Way¬ 
like galaxies. The circular velocities of non-relativistic 
objects are smaller than the GR prediction at fixed ra¬ 
dius and light is less strongly lensed. We discuss our 
results and conclude in section jV] The metric signature 
is (—, +, -b, -b) and the Planck mass Mpi^ = (SttGn)”^. 


II. MODIFIED GRAVITY THEORIES AND THE 
VAINSHTEIN MECHANISM 


Theories that screen using the Vainshtein mechanism 
are abundant in the literature (see | 3 ] for example) and 
galileons are the most common paradigm by far. These 
are defined in the decoupling limit by the following La- 
grangiar 0 
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Here .. .V and X = -l/2d^4>d>^(^ 

The theory we will study below does not have a confor¬ 
mal coupling a but instead sources the galileon through 
couplings to the curvature tensor present when the quar¬ 
tic term is covariantised. 

Despite being higher-order in the action, the equations 
of motion are precisely second-order and so there are no 
ghost instabilities. Whilst fine around flat space, the final 
two Lagrangians do not give rise to second-order equa¬ 
tions of motion about an arbitrary space-time because 


^ We ignore the tadpole term here. 

® Note that our definition of X differs from 14211431 by a factor 
of —2 but conforms with m- Furthermore, our fields are di¬ 
mensionless, which conforms with neither previous work but has 
the advantage that the equation of motion for 0 is similar to the 
Poisson equation. 
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their variation leads to couplings of the field to curvature 
tensors that are identically zero when the space-time is 
MinkowskQ These terms always contain third-order (or 
higher) derivatives of the metric or fields. In fact, the 
Lagrangians above are not unique and there is an infi¬ 
nite number of equivalent descriptions that are related 
via integration by parts. Indeed, one can remove the fi¬ 
nal term in £3 and the terms in the curly brackets in 
£4 and £5. One way to restore the ghost-free nature 
of the theory is to start from these simpler forms and 
add additional couplings of the fields to curvature ten¬ 
sors so that the higher-order derivatives exactly cancel. 
With this procedure, one is led to covariant galileons [HH] . 
These are a subset of the Horndeski Lagrangian [41] , the 
most general Lagrangian with second-order equations of 
motion. Recently, jH] |33| (see also [H]) have shown 
that there exist Lagrangians that give rise to higher-order 
equations that are still ghost-free due to the existence of 
hidden constraints. The naive covariantisation of £4 and 
£5 are a subset of the most general Lagrangian and so 
this presents us with a new and interesting theory that 
screens using the Vainshtein mechanism that has yet to 
be explored. In fact, one can show that the static, spher¬ 
ically symmetric Vainshtein solution is unstable to per¬ 
turbations when C 5 ^ 0 [bOj and so from here on we will 
investigate the Lagrangian 


£ 

7^ 


= Mpi2 


R £4 

2+^+AJ 


( 6 ) 


P5] have dubbed the entire class of beyond Horndeski 
theories and so we will refer to this theory as the 
G^-galileon. 

Before examining the breaking of the Vainshtein mech¬ 
anism in this covariantisation, we briefly review the 
screening of the quartic galileon in the decoupling limit 
for comparative purposes. The equation of motion 
for static, spherically symmetric over-densities, = 
diag(p(r), 0 , 0 , 0 ), is 


j__d 

r2 


+ A4' 


j,/3 


SairGp. 


(7) 


The first term is simply and is due to the quadratic 
term in the Lagrangian, whilst the second term is due 
to the quartic Lagrangian. The fifth-force £5 due to the 
conformal coupling (see gSlIlTI for a discussion of the 
disformal coupling) is £5 = —a(j)' and so one can inte¬ 
grate Q to find the algebraic relation 


q,2^4^2 


Fi = 2a" £n, 


( 8 ) 


where £n is the Newtonian force. One can see that on 
large enough length scales (to be made precise shortly) 


^ This is due to the relation required 

when varying with respect to the field. 


£5 = 2 a" £n is a consistent solution of the equation and 
the force is unscreened. On smaller length scales, one has 


Fn 



r-v = 2 a 


GM 


(9) 


where we have assumed a point mass source p{r) = 
M(5(3)(r) in order to provide a concrete expression. The 
Vainshtein radius ry defines the transition between the 
screened and unscreened regime. When r ^ ry, the 
force is unscreened and when the converse is true one 
has £5 ^ £n and the force is screened. Typically, an 
object characterised by some length scale R has a Vain¬ 
shtein radius ry ^ R, for example, the Vainshtein radius 
of the Sun encompasses the entire solar system [5ll|52|. 
For this reason, the Vainshtein mechanism is very effi¬ 
cient at screening and there have been few astrophysical 
signatures to date (see [22l |53l |54| for some exceptions to 
this). Below, we will examine the non-relativistic limit of 
the G^-galileon and show that the Vainshtein mechanism 
is only partially effective. 


A. The Non-Relativistic Limit 


We wish to derive the equations governing Newto¬ 
nian perturbations about an FRW space-time sourced 
by some non-relativistic density profile p{r). To this 
end, we expand the field about its cosmological value 
4>(r, t) = tpoif) + ’’■(j'j i) and the metric as 


ds" = — [I -I- 2$(r, t)] dt^+a‘^{t) [I — 24'(r, t)] dij dx* da;". 

( 10 ) 

The radial coordinate is r = ay/6ijX'^xF In what follows, 
we will assume that the local effects are due primarily to 
the cosmological boundary condition on </) and not due 
to effects coming from the FRW metric. To this end 
we will ignore all factors of a{t) and Hitf\ We also 
ignore any time-dependence since we are interested in 
static situations in this work. Next, we must solve the 
equations of motion for the dynamics of $, 4' and tt. 

This has already been studied by [IS] for the full set 
of beyond Horndeski theories. There, they expanded the 
action to the relevant order in metric and field perturba¬ 
tions in order to derive an effective non-relativistic action. 
One must be careful when gauge-fixing the action. There 
are certain criteria for doing so in order that information 
is not lost in the resulting equations of motion and, as 
discussed by , fixing the conformal Newtonian gauge 
(101 at the level of the action may result in equations 


® This is an ad hoc assumption that must be tested for self- 
consistency. In fact, if one were to include these factors, the 
only change would be in how the final new parameter T is de¬ 
fined in terms of the parameters appearing in the fundamental 
theory. Since we will always work with T and not bo or A we are 
not missing any potential new effects by making this assumption. 
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of motion that are not equivalent to those found from 
the gauge fixed forms of the equations found using the 
full action. For this reason, we have derived the non- 
relativistic limit by gauge-fixing the equations of motion 
and not the action. This derivation can be found in Ap¬ 
pendix In fact, our equations are identical to those 
of |45j when their parameters are chosen to match ours 
but there is no underlying reason for this and, indeed, it 
is not necessarily the case that this is the case for more 
general models. 

In what follows, we will use the variables (which are 
similar to those used by |i5] l 


X = 




y = 


r 


M{r) 

SnMpi^r^ 


( 11 ) 


where a prime denotes a derivative with respect to the 
radial coordinate and M{r) is the mass enclosed within 
a sphere of radius r. Using these variables, one finds the 
following equations in the non-relativistic limit (after in¬ 
tegrating the equations of motion by parts, see Appendix 
0: 


z-Hya;2 = A, (12) 

y — z — — ex [rx' -|- x) = 0 and (13) 

x^ 

X — 4^ -I- lOexy + 2ex (2z + rz') = 0, (14) 


where 


e = 



(15) 


One can then use equations (12) and (13) in equation 
(14) to find an algebraic equation for x: 


second-order. Taking A ^ ij^one find^ 

4 7A -|- r A' 


= £A^ 


10T-b2 


At this point we introduce the new parameter 


T = = 



(17) 


(18) 


This is the parameter that quantifies the size of the new 
effects presented here and we will work with it exclusively 
from here on. Using (17), the solutions for y and z are 


y = 


z = 


(r^A)' 


l-b5T 4(1-b5T) r 


and 


A 


5T 




l-b5T 4(1-b5T) r2 ' 

One would ideally identify Newton’s constant as 

G 


Gn = 


1-f 5T’ 


(19) 

( 20 ) 

( 21 ) 


but this definition is subtle since we do not know the 
solution for the metric potentials outside of the source. 
Traditionally, one works in a gauge where, to Newtonian 
order. 


= - (l + + (l - 2,^5®) ,22) 


outside the source, which defines Gn- We will show in 
Appendix [b] that (21) is indeed the correct definition of 
Gn- This then allows us to write (19) and (20) in the 
more useful form 


d$ 

dr 

GNM(r) , , 

= ,.2 +4 GnM (r) 

(23) 

d4' 

GNM(r) 

5T GNM'(r) 

(24) 

dr 

j,2 

4 r 


These will be our starting point when deriving the new 
properties of stars and galaxies- Before doing so however, 
it is worth pausing to note three important properties of 
these equations- First, the extra terms that appear are 
not due to contributions from x ~ tt'; these terms have 


-f 20s^^ x^ - 2e [(7A -b rA')] -1 = 0- (16) 

Note that x' does not appear in this equation; this is the 
hidden constraint of the beyond Horndeski action. De¬ 
spite the fact that the equations of motion appear third- 
order, only two integrations are required to fully solve 
the system; the hidden constraint has allowed us to find 
a suitable linear combination of the equations that are 


® We assume rA' ~ A. _ 

^ The careful reader will notice that l |l6| is in fact a cubic equation 
before cancelling a common factor of x and that x = 0 appears 
to be a consistent solution. In fact, there is a small constant 
term in the cubic equation that is Hubble suppressed m and 
so this solution cannot be realised in practice. Furthermore, 
we are using the G^-galileon as an example of the breaking of 
the Vainshtein mechanism in beyond Horndeski theories. The 
equations we will ultimately obtain are correct for any beyond 
Horndeski theory and are obtained by taking the x ^ 0 branch 
of solutions. See m for the general derivation. 
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been suppressed by the Vainshtein mechanism. Instead, 
they are due to the cosmological time-derivatives that the 
Vainshtein mechanism fails to suppress. This behaviour 
has previously been noted in the case of covariant quartic 
galileons [ID], however there the effects were simply a 
time-dependent rescaling of Newton’s constant. Second, 
we will show in Appendix that the PPN parameter 
7 = 1 and this means that, a priori, there are currently no 
constraints on the parameter Tfj^ In practice, we expect 
it to be small and one of the aims of this work is to discern 
the typical values which lead to large deviations from GR. 
Third, the final term in equation (23) is typically negative 
provided that the density decreases outwards (as it does 
in most astrophysical systems). Since d$/ dr governs the 
motion of non-relativistic particles, the theory predicts 
that any such objects moving inside a larger source will 
feel a reduced strength of gravity than that predicted by 
GR. 


III. STELLAR STRUCTURE TESTS 


In this section we derive the new stellar structure 
equations predicted by the theory and use them to pre¬ 
dict novel phenomena that may be used as observational 
probes. We note that, while the derivation of the equa¬ 
tion for dd)/ dr used the G^-galileon as its starting point, 
the final form applies to any beyond Horndeksi theory 
pp] . All that changes is the definition of the parameter 
T in terms of the fundamental parameters that appear in 
the Lagrangian. For this reason, equation (23), and the 
results we obtain here, are completely general and one 
may treat T as a free parameteiQ 

Since we are working in the Jordan frame, the energy- 
momentum tensor of matter is covariantly conserved and 
so one has = 0. Taking the energy-momentum 

tensor to be that of a static, spherically symmetric non- 
relativistic source, the Euler equation found from this 
conservation law is 


1 dP 
p dr 


d$ 

dr 


(25) 


using equation (23), one finds the new hydrostatic equi¬ 


librium equation, which governs the structure of stars: 


^ = (26) 


® It is known that the theories here are related to Horndeski the¬ 
ories via a disformal transformation m and so the constraints 
of 15611571 may apply to T once one has accounted for the Vain¬ 
shtein mechanism and specified a cosmology. 

® This is because non-relativistic objects such as main-sequence 
stars respond to <I> only. In the most general theory, the param¬ 
eter appearing in the equation governing 'll may differ from the 
one appearing in the equation governing 'I', and so one can have 
two free parameters in general. 


This is the only modification to the stellar structure equa¬ 
tions in this theory. The conservation of mass is unal¬ 
tered: 


= 4:Trr‘^p{r) (27) 

dr 

and for later convenience one can differentiate this to find 


d^M 

dr 2 


Snrp + 47rr^ 


dr 


(28) 


The other two stellar structure equations do not depend 
on the theory of gravity; they are the radiative transfer 
and energy generation equations respectively [58] : 


dr 3 K pL 

dr 4a Airr'^ 

= 47rr^e(r), 
dr 


(29) 

(30) 


where L is the luminosity, k is the opacity, e(r) is the en¬ 
ergy generation rate per unit mass and a is the constant 
relating the radiation pressure to the temperature (see 
below). Equation (26) is all that is required to study 
the structure of non-relativistic stars in these theories. 
We will do this below, first for simple polytropic mod¬ 
els and then for realistic stars found using a consistent 
numerical implementation of the new hydrostatic equi¬ 
librium equation. The former models are not powerful 
enough to allow a comparison with real data but they 
are useful in that they decouple the non-gravitational 
physics, which allows one to gain an analytic handle on 
the new effects arising in these theories without the com¬ 
plications arising from thermodynamic effects. Further¬ 
more, main-sequence stars are well-approximated by cer¬ 
tain Lane-Emden models and this will allow us to discern 
the new features one would expect in stars like the Sun. 
Lane-Emden models have been particularly useful in un¬ 
derstanding how stars behave in chameleon (and other 
similar) theories of gravity [TDl [TD| |TD| . 


A. Polytropic Models 

In the absence of any knowledge of the micro-physics, 
the stellar structure equations do not close and one must 
supply and equation of state relating P and p. Polytropic 
models are defined by the equation of state 

P = Kp'^. (31) 

n is known as the polytropic index. Main-sequence stars 
are well-described by n = 3 and fully convective stars 
are well-modelled by n = 1.5. In GR, polytropic models 
lead to self-similar equations of stellar structure and since 
the new parameter T is dimensionless the same is true 
in this theory. We can then scale out the dimensionful 
quantities in the theory by writing r = Tc^ p = /Oc0(C)" 
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and P = , where Pc and pc are the central 

pressures and densities and 


(n + 1)P c 
AttGnPc'^ 


(32) 


Note that Pc and pc are related via Pc = . 

One can then insert equation (31) into equation (26) and 


use equations (27) and (28) to find 


1 d 


d0 T 




= -r. (33) 


When T = 0 this equation describes the structure of 
polytropic spheres collapsing under gravity described by 
GR. In this case, it is known as the Lane-Emden equa¬ 
tion and so we refer to ( [3^ and the modified Lane-Emden 
equation (MLE). The boundary conditions for the dimen¬ 
sionless density 9 are 0(0) = 1 (so that p(0) = pc) and 
0'(O) = 0, which is a requirement of spherical symmetry. 
The radius of the star is found from the value of de¬ 
fined by 9{^r) = 0 so that R = Vc ^ r - One can then see 
that the extra terms in equation (33) vanish at the sur¬ 


face of the star and for this reason, stars in this theory 
share many features with GR. In particular, the mass of 
the star is 

cR 

M= drAirr^ p{r) 

Jo 

r^R 


47rrc^Pc [ = AttTc^PcUJr, 

Jo 


where 


W_R = -^R 


d9 


(34) 


(35) 


This is precisely the same formula that is found in GR. 
The only difference being that ujr will be different de¬ 
pending on the value of T. In particular, this means that 
whereas a star with fixed central density will have dif¬ 
ferent masses between the two theories, the mass-radius 
relationship R oc predicted by GR is unal¬ 

tered. 


One can in principle solve equation (33) to find the new 
properties of many polytropes. The goal of this paper 
is to predict the new observational features this theory 
predicts and so we will now specialise to the case of main- 
sequence stars. 


1. Main Sequence Stars 

We will investigate the new properties of these stars 
using the Eddington standard model. The key assump¬ 
tion of this model is that the specific entropy gradient 
(a T^/p) is constant, which allows us to decouple the 
radiative transfer and energy generation equations. The 
opacity is taken to be constant and we assume that the 


pressure supporting the star is due to the thermodynamic 
properties of the gas. 


p — 

Rgas — 


pk-gT 
pmn ’ 


(36) 


where p is the mean molecular mass, and the radiation 
pressure from photons generated by hydrogen burning in 
the core, 


Erad = 


so that the total pressure is 


pma 3 


(37) 


(38) 


Defining the quantity fi = Pgas/P, one can equate Pgas 
with Prad/(1 — fi) using (36) and (37) to find 


1-/3 Sknp 




(39) 


The assumption that the specific entropy gradient is con¬ 
stant ensures that ,5 is a constant, in which case one can 
use equation (39) in equation (1M| to find 


P = K{P)p^ 


with 


K{P) = 


k-Q 

pmu 


)'(' 


-P 


(40) 


(41) 


Equation (40) is the equation of state of an n = 3 poly¬ 


trope and so we can describe these stars by solving the 
modified Lane-Emden equation. One important observa¬ 
tional property of main-sequence stars is their Luminos¬ 
ity, which we can predict as follows. Differentiating equa¬ 
tion (371 and substituting it into the hydrostatic equilib¬ 


rium equation (setting P = (1 —/3)Prad) one finds (at the 
surface of the star) 


L = 


47rGNM(l - P) 


(42) 


This is identical to the expression found in general rel¬ 
ativity and so any observational differences (at fixed 
mass), if any, are encoded in p. One then requires an 
expression for p. This can be found as follows: Inserting 
the definition of Vc (p^ and equation (41) into equation 


(34) one finds the modified form of Ec 


dington’s quartic 


equation: 


where 


1-/3 

/34 


(gfiJii'i " 


M'Edd = 


dujR f k-B 






Hi 


(44) 
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is the Eddington mass and Cor « 2.018 is the value of 
lor when T = 0 (i.e. the result in GR). Numerically, one 
has MEdd ~ 18.3^“^ and from here on we take = 0.5, 
corresponding to fully ionised hydrogen, which is appro¬ 
priate for main sequence stars burning on the PP-chains. 
Equation (431 then shows that all of the modified grav¬ 
ity effects are encoded in the new value of ujr, which is 
a number describing the change in the structure of the 
star. 

We are now in a position to integrate the MLE nu¬ 
merically and derive the new properties of main-sequence 
stars. It is well-known that there are subtleties with inte¬ 
grating the Lane-Emden equation numerically and these 
persist with the MLE. We briefly describe how the equa¬ 
tion is integrated numerically in Appendix We begin 
by plotting the Lane-Emden function 0{^) for various val¬ 
ues of T ~ 10“^ in figure The deviations from the GR 
curve are indistinguishable when T < 10“^. One can 
see visible changes in the shape of the function and the 
intersect on the ^-axis (i.e. the value of ^r), namely that 
solutions with larger T have larger values of Indeed, 
plotting uir/lor and ^r/^r as a function of T in figure 

shows that one expects large signals when T > 10“^. 
That being said, one should not draw too many conclu¬ 
sions directly from 9{^) and ^r. 



FIG. 1. The n = 3 solutions of the Lane-Emden equation 
(black) and the MLE with T = 0.1 (blue) and T = 0.2 (red). 


Using these solutions, we can solve equation ( |43[ ) to 
find /3(M), which can be used in equation 


to find 

the ratio of the luminosity predicted in the new theory 
to that predicted by GR at fixed mass. In figurej^we plot 
this ratio for a IMq star as a function of T. Once can 
see that there is a visible decrease in this ratio around 
T ~ 10“^. In figure]^ we plot the same ratio as a func¬ 
tion of mass for three fixed values of T. One can see 


One should be careful to note that this does not necessarily im¬ 
ply that the stars have larger radii since this is degenerate with 
changing the central density and K. Indeed, if one wishes to 
compare two stars then it is important to decide which quan¬ 
tities to fix and which to vary in order to make a meaningful 
comparison. 



FIG. 2. ljJrIujii (red) and ^r/^r (blue) as a function of T for 
n = 3 solutions of the MLE. 


that the modifications affect low mass stars to a greater 
extent than high mass stars. This is because low mass 
stars are gas pressure-supported whereas high mass stars 
are radiation pressure-supported. The luminosity of gas 
pressure-supported stars is more sensitive to the strength 
of gravity than radiation pressure-supported stars (see 
m for a discussion on this), hence the turn towards GR 
at high masses. The stars in the new theory are dim¬ 
mer by a factor of 5 for T ^ 0.2 and the ratio increases 
steadily towards 1 when T is decreased. In particular, the 
ratio is greater than 0.9 when T < 0.1. Very dim stars 
of IMq are clearly at odds with observations of the solar 
brightness and so one can tentatively place the constraint 
T < 0(10- Said another way, since the cosmology 
of this theory has yet to be worked out, any theory with 
a best-fitting model that requires T > (9(10“^) is likely 
in tension with astrophysical constraints. 


This is only tentative because polytropic models are subject to 
the caveats mentioned above and these effects may be diminished 
by or be degenerate with non-gravitational physics not included 
in our model. Furthermore, it may possible to achieve the ob¬ 
served solar luminosity in more complicated models by changing 
unmeasured quantities such as the solar equation of state or the 
neutrino loss rate. Indeed, we will see in the next section that 
the effect of increasing T is degenerate with increasing the metal- 
licity. 













^MG 

t-GR 



FIG. 3. The luminosity enhancement for a IMq star as a 
function of T. 
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FIG. 4. The luminosity enhancement as a function of mass 
for T = 0.3 (red), T = 0.2 (blue) and T = 0.1 (green). 


B. Realistic Stars 

The simple semi-analytic model above suffers from sev¬ 
eral drawbacks. Firstly, it does not include any non- 
gravitational physics. Secondly, it predicts that the sur¬ 
face temperature, one of the most important observa¬ 
tional properties of stars, is zercp] and so the modifica¬ 
tion to the stars colour, if any, is unknown. Finally, it 
is unable to treat post-main-sequence stars such as red 
giants due to their complex composition and the absence 
of any convective physics in the model. 

For these reasons, we have modihed the publicly avail¬ 
able code MESA [3S1 SS] to solve the modified stellar 


One can see from equation l |39| that T oc 6 and since the star’s 
radius is defined by 0(^r) = 0 the temperature is identically 
zero at the stellar radius. More realistic stellar models define 
the effective temperature as the temperature of the photosphere, 
which is the radius at which the optical depth is 2/3. 


structure equation (261. MESA is a self-consistent stellar 
evolution code that solves the stellar structure equations 
coupled to convection recipes, nuclear reaction networks, 
atmosphere models and radiative transfer processes. It 
is able to produce stellar models that are realistic and 
can be compared with observational data. Indeed, it 
has already been used to place the most stringent con¬ 
straints on chameleon theories [iniiiiiis] to date. Grav¬ 
itational physics only appears in the MESA code via the 
hydrostatic equilibrium equation. All other equations en¬ 
code the effects of processes such as convection or nuclear 
burning, which are not directly controlled by the theory 
of gravity. For this reason, the only change we have made 
is to replace the GR equation with equation (26). 

As an example, in figure we plot the Hertzprung- 
Russell (HR) diagram for a solar mass star. Included are 
the tracks for a star with solar metallicity {Z = 0.02) as 
predicted by GR and G^-galileons with various values of 
T indicated in the caption. One can see that the main 
effect is in the temperature rather than the luminosity. 
Main-sequence stars have a smaller effective temperature 
than stars at the same point on the HR track when the 
modification is stronger. An examination of the tracks 
shows that the effect is far more pronounced in main- 
sequence rather than red giant stars. As with any as- 
trophysical system, one must worry about degeneracies. 
Stars with higher metallicities have smaller effective tem¬ 
perature^^ at fixed luminosity and this is one potential 
source of degeneracy. Indeed, figure also includes the 
GR tracks for stars with Z = 0.01 and Z = 0.03. One 
can see that the main-sequence track for a G^ star with 
Z = 0.02, T = 0.1 is nearly coincident with the main- 
sequence track for a GR star with Z = 0.03. On the 
main-sequence, the effect of G^-galileons is to mimic stars 
of larger metallicity. This is not the case for post-main 
sequence stars because the tracks converge to a value very 
close to the GR prediction whilst the higher metallicity 
GR star still has a smaller effective temperature. Finally, 
we indicate the point on each track when the age of the 
star is solar age. One can see that stars with stronger 
modifications evolve more slowly and are less ephemeral 
than their GR counterparts (at fixed metallicity). 

As a final example, we plot the HR tracks for 2Mq 
stars in figure!^ Shown are the tracks for Z = 0.02 stars 
in GR and G^alileons with T = 0.3 and T = 0.1. One 
can see a similar effect to that found for solar mass stars 
where the effective temperature is smaller with increas¬ 
ing T but the decrease in the luminosity is also more 
pronounced. This is in contrast to what we found in the 
previous section. Higher mass stars have a composition 
that is far less uniform than low mass stars and are not 
as well-described by polytropic models. If one is simply 
changing the value of Gn then the change in the lumi¬ 
nosity is due entirely to the change in /3, however, since 


This is because they have more free electrons and the opacity is 
larger. 
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FIG. 5. The HR diagram for a solar mass star. The red curve 
is the GR track when the metallicity is solar {Z = 0.02) and 
the blue and green tracks show the evolution of the same star 
when T = 0.3 and T = 0.1 respectively. The black dashed 
(upper) and dotted (lower) lines show the tracks for the same 
star when the theory of gravity is GR and Z — 0.01 and 0.03 
respectively. The black circles indicate the point on the track 
when the age of the star is 4.6 x 10® yr. This is not shown for 
the Z — 0.01 and 0.03 stars. 


equation (26) depends primarily on the internal struc¬ 
ture of the star it is possible that polytropic models miss 
important contributions coming from non-uniform com¬ 
positions. Finally, we have also plotted the track for a 
GR star with Z = 0.03. One can again see that the main- 
sequence portion of this track is nearly coincident with 
the G^-galileon star of solar metallicity, showing that the 
degeneracy is still present for higher mass stars. Again, 
the two tracks differ in the post-Main-sequence phase. 



FIG. 6. The HR diagram for a 2Mq star. The red curve is 
the GR track when the metallicity is solar {Z — 0.02) and 
the blue and green tracks show the evolution of the same star 
when T = 0.3 and T = 0.1 respectively. The black dashed 
line shows the track for the same star when the theory of 
gravity is GR and Z = 0.03. 


profile. In what follows, we will investigate two potential 
signatures of the breaking of the Vainshtein mechanism: 
rotation curves and lensing. In all cases we will assume 
that the dark matter is distributed according to the the 
Navarro-Frenk-White (NFW) profile |59) : 



We will always examine Milky Way-like galaxies and so 
we will take rg = 20 kPc and ps = 6.68 M 0 kPc~^ in order 
to provide concrete numbers. These are typical values 
for the Milky Way (assuming a halo mass of 10^^Mq) 

mM- 


A. Rotation Curves 


In general relativity, the rotation curve of a galaxy is a 
direct probe of the derivative of the gravitational poten¬ 
tial and is hence a probe of the interior mass distribution. 
Here, we will take the density profile as giveij^ and use 
the modified equations for $ and dl to predict the new 
galactic properties given this profile. 

Assuming circular motion, the velocity and radius of a 
test body inside the halo are related via 


_ d$ 

r dr 


and so using equation (231 one finds 


(46) 


47rGN?^s^Ps 

r 



4 (1 + ^)3 
(47) 


The rotation curve for a Milky Way-like galaxy is plot¬ 
ted in figure for various values of T indicated in the 
caption. One can see that the discrepancy between the 
two theories is minimal at small radii and so the part of 
the curve measured using stellar velocities is largely un¬ 
changed. The discrepancy is more pronounced at larger 
radii where the curves begin to flatten. Here, G^ the¬ 
ories with larger values of T predict that objects orbit 
with smaller circular velocities than predicted with GR. 
The velocities in this portion of the curve (in the Milky 
Way) are typically measured using HI emission lines. 


B. Strong Lensing 


The bending of light inside a source is governed by the 
lensing potential (4>-|-4'). Using equations (23) and (24) 


IV. GALAXY TESTS 


Galaxies are another astrophysical object that can po¬ 
tentially exhibit effects due to their extended dark matter 
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One should ideally run new simulations in this model to verify 
that NFW is still a good fit but this is clearly beyond the scope 
of the present work. 
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FIG. 7. The rotation curves for the NFW profile with Tb = 20 
kPc and ps = 6.68. The GR prediction is shown in black and 
the red and blue curves correspond to T = 0.3 and T = 0.5 
respectively. 
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FIG. 8. The lensing potential predicted by the NFW profile 
with ra = 20 kPc and ps = 6.68. The GR prediction for 
this ratio is unity. The G^-galileon predictions with T = 0.1, 
T = 0.3 and T = 0.5 are plotted using black, red and blue 
curves respectively. 


one finds 


) 

2 ■ 

(48) 

In figurej^we plot the ratio ($ + 'I')/2<I>, which quantifies 
the amount by which light is bent by gravity relative 
to the gravitational force felt by objects moving at non- 
relativistic velocities, for various values of T indicated 
in the caption. GR predicts that 41 = $ so that this 
ratio is unity. One can see from the figure that this ratio 
is less than unity in G^-galileon theories and decreases 
with increasing values of T. This means that the lensing 
mass is smaller than the dynamical mass and so strong 
lensing could be used as a novel probe. Indeed, |62] have 
already used strong lensing, which probes the weak field 
metric on kilo-parsec scales, to constrain deviations of 
the parameter 7 , defined by 


$ + \]/ = - 


SttCnTs^Ps 


In 


rs 


(1 


ds^ 



dt^ + 


1 - 27 



da;^ (49) 


from unity to the 5% level. Their constraint cannot be 
used to constrain our model parameters since <f> and 
are not inverse-square in our model. Nonetheless, this 
serves as an example of how astrophysical measurements 
such as strong lensing can probe this theory of gravity. 


V. DISCUSSION AND CONCLUSIONS 

In this paper we have studied the non-relativistic limit 
of a subset of beyond Horndeski theories, the covarianti- 
sation of the quartic galileon or, as we have dubbed it, the 
G^-galileon. We have found that whilst the Vainshtein 
mechanism is active in suppressing gradients of the scalar 
degree of freedom, the equations of motion for the met¬ 
ric potentials are different from the Poisson equation. In 


particular, there are extra terms inside sources of finite 
extent proportional to the first and second derivatives 
of the mass. These terms act to weaken the strength of 
gravity provided that the density decreases outward, as 
is the case for most astrophysical sources. 

This opens up the possibility of testing these theories 
using the astrophysics of finite sources and in this pa¬ 
per we have taken the first steps towards identifying new 
and novel probes. The strength of the deviations can be 
quantified by a single number, T, which is a combina¬ 
tion of the time-derivative of the cosmological scalar and 
the new mass scale A appearing in the quartic galileon 
Lagrangian. To-date, the specific theory has yet to be 
constrained using other, more conventional probes, and 
so, a priori, there are no constraints on T. 

One promising probe is the equilibrium structure of 
stars. Indeed, these have already been used to place the 
most stringent constraints on chameleon-like theories to- 
date. We have derived the modified hydrostatic equi¬ 
librium equation and have used it to study both simple 
semi-analytic polytropic stars and fully realistic stars by 
updating the stellar structure code MESA. In the former 
case, the modified Lane-Emden equation was derived for 
a general polytropic index before specialising to the case 
of main-sequence stars. The weakening of gravity causes 
stars to be dimmer than their GR equivalents at fixed 
mass. Full MESA simulations revealed that the surface 
temperature is also lower and so one would generally ex¬ 
pect stars of fixed mass to be dimmer and redder than 
their GR counterparts. We also found that the G^-stars 
evolve at a slower rate than the equivalent GR stars. As 
with any astrophysical test of gravity, there are degenera¬ 
cies with non-gravitational astrophysical processes. We 
identified a degeneracy with metallicity where a G^-star 
of fixed mass and metallicity mimics a GR star of the 
same mass with a larger metallicity on the main sequence. 
This is not true on the red giant branch, where the G^- 
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star tracks are nearly identical to the GR tracks for stars 
of the same metallicity. This then suggests a potential 
probe of these theories whereby the metallicity (and pos¬ 
sibly age) of a globular cluster inferred from the main- 
sequence portion of the HR diagram would differ from 
that inferred using the red giant branch. 

Next, we turned our attention to galactic probes. We 
calculated the rotation curves for Milky Way-like NFW 
halos and found that non-relativistic objects orbit with 
smaller circular velocities in G^-galileon theories com¬ 
pared with GR. The effect is minimal on the portion of 
the curve that rises steeply but can be large when the 
curves flatten. The leasing potentials were also calcu¬ 
lated and it was found that the leasing mass is smaller 
than the dynamical mass. The difference can be at the 
10 % level or greater for certain parameter choices and so 
this could potentially be probed using strong leasing. 

Before concluding, it is worth commenting that the 
specific model studied here is known to have a gradient 
instability [53] at late times. Whereas this is not catas¬ 
trophic for the theory, it is worrying. One can cure this by 
adding more Horndeski (and possibly beyond Horndeski) 
terms. In terms of our calculations here, this changes 
very little. Following the analysis of [IS], no new terms 
in the equation of motion for $ are generated and all that 
changes is the definition of T in terms of the fundamental 
model parameters. Our predictions for non-relativistic 
objects (stars and rotation curves) are therefore robust 
to this change since we have quoted all of our results in 
terms of T. The only difference is the introduction of a 
new parameter in the equation of r noti on for 4', by which 
we mean T —> T 7 ^ T in equation (24) whilst T remains 
in (23). In this case, the leasing potential would depend 
on two parameters but would look qualitatively the same. 
This is only the case when the Lagrangian contains quin- 
tic terms (either Horndeski or beyond Horndeski) and 
couplings of the scalar to the Ricci scalar proportional to 
X (see US], Appendix 2). 

The breaking of the Vainshtein mechanism is only 
present when the cosmological scalar has a non-vanishing 
time-derivative. All of the effects presented here were 
non-negligible when T > 0{1) and so any theory whose 
best-fit cosmology predicts values of T in this range can 
potentially be probed using astrophysics. 


Finally, we conclude by discussing future avenues for 
testing the breaking of the Vainshtein mechanism. Here 
we have only studied the equilibrium structure of stars 
but this is only probed using the surface properties. Stel¬ 
lar oscillations and helioseismology allows one to probe 
the structure of stars and contains far more informa¬ 
tion. Indeed, the strongest constraints on f{R) theories 
comes from stellar pulsations m- Since the theory pre¬ 
dicts that our own Sun is subject to the breaking of the 
Vainshtein mechanism, one potential probe would be to 
look at radial and non-radial modes of stellar oscillations. 
We have also seen that the tracks in the HR diagram 
for G^-stars mimic those in GR with higher metallicities 
and so one potential signature may arise if larger mass 
stars cross the instability strip at different locations at 
fixed metallicity. This would manifest as a change in the 
period-luminosity-mass-colour relation for Gepheid and 
RR-Lyrae variable stars. Another interesting possibil¬ 
ity is convective stars, which have n = 1.5. Examining 
the modified Lane-Emden equation reveals that, for these 
stars, it contains a term proportional to which 

diverges at the stellar radius. This may suggest that fully 
convective stars may not exist in these theories. White 
dwarf stars are also subject to the modifications and it 
is possible to extend our polytropic analysij^ to derive 
their new properties and the Chandrasekhar mass to see 
if it is modified. 

All of these effects could potentially probe the entire 
class of beyond Horndeski models and here we have pre¬ 
sented the tools required to do so, using the G^-galileon 
as a specific example. 
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Appendix A: Derivation of the Non-relativistic Limit 

In this Appendix we derive the equations of motion 
for the scalar field perturbation tt and the gravitation 
perturbations $ and 'I' in the non-relativistic limit fol¬ 
lowing the approach developed in the Horndeski theory 
[501164] . We first ignore terms involving time-derivatives 
of these fields. We then assume that the perturbations 
are small and neglect non-linear interactions containing 


higher-order powers of these fields and their first-order 
derivatives. On the other hand we keep all terms with 
second order derivatives that preserve the Galileon sym¬ 
metry. In addition, unlike the Horndeski theory, the 
equation of motion for tt in our theory involves a third 
order derivative of 4/, which breaks the Galileon symme¬ 
try. The scalar field equation and the Einstein equations 
give 


^ (V^7r)2 - (V*Vj7r)(V*V^'7r) 


= 47rGp, 


- 


V^TT- 


41 
2 ■ 
A4 . 


(V^tt)" -b 3(V,V^-7r)(V*V^7r) -b 4(V,V^7r)(VV) 


= 0 , 


(VV)^ - 3(V^7r)(V,Vj7r)(V*V^7r) -b 2(V,Vj7r)(VfcVV)(V'=W) 


(Al) 

(A2) 

(A3) 


£ 5(V2$)(V27r) - 5(ViV^-$)(V*V%) -b + (V,Vj4')(V*V%) -b 2(ViV2^')(W) 


= 0 . 


Imposing the spherical symmetry and integrating by 
parts, these equations reduce to equations (12), (13) and 

pl. 


Appendix B: The Vainshtein Mechanism Outside 
Finite Sources 


In this Appendix we will examine the exterior profile 

of the metric potentials 13 

We begin by discussing the boundary conditions. The 
assumptions of spherical symmetry require us to impose 
the boundary conditions 7r'(0) = $'(0) = 4''(0) = 0. 
Equations (14)-(13) were found by integrating a set of 
second-order equations once and dividing by r^. In the¬ 
ory, these contain arbitrary integration constants on the 
right hand side of the form Cijr^, but since these diverge 
as r —0 we must have Gi = 0 and so these terms are 
absent. The same is not true outside of the source. The 
exterior solution can contain terms that scale as Bijr^ 
since the solution is valid for r > R only. The boundary 
conditions at the stellar radius fix these integration con¬ 
stants. Now R is typically defined such that p{R) = 0. In 
section we found that we had to impose two bound¬ 
ary conditions at the centre of the staip^ and since the 
equations are second-order, this leaves no freedom to fix 
any conditions at the stellar surface. In particular, we 


We would like to thank Eugeny Babichev, David Langlois and 
Ryo Saito for useful discussions, which helped us identify the 
incorrect use of the boundary conditions an earlier verison of 
this work. 

This is true in general. One requires dP/ dr(0) = 0 due to spher¬ 
ical symmetry and the central pressure must also be specified. 


are not free to set p'{R) = 0 by hand. We will use the 
notation p' = p'[R). 

Accounting for the three integration constants Bi, the 
equations of motion outside the source are 


5 £ 2 Bi 

Z + —X - T = 0, 

2 r® 

£ 2 / / \ B 2 

y — z - X — ex irx -b x)-= 0 

2 r® 


— 4^ -b IQexy -b 2 £x (2z -b rz') — ^ = 0. 

A4 j-i 


(Bl) 

and 

(B2) 

(B3) 


After some algebra, one can combine these to find a single 
equation for x: 

x-4(I-b 5T) ^-b(8Hi-b IOB 2 ) ^^-b^ = 0. (B4) 

This is a cubic equation for x, which in general has three 
distinct solutions. 


Since equations (B1)-(B3) do not contain any 15- 


function sources we expect that B^ = 0 both inside and 
outside the source to avoid discontinuities in x and x'. 
We therefore make the ansatz that B^ = 0 and check 
that this indeed leads to the correct solution that satis¬ 
fies all of the boundary conditions. Thus, the solution is 
self-consistent. Solving equation (B4) gives 


x2 = 


4Bi -b 5 B 2 

' 2 -b lot ' 


^^+G(r3/2). (25) 


Equations (iBlj) and (B2) contain x' and z' and so we 


must require that x and z are continuous at the stellar 
radius. This fixes Bi and B 2 as 


Bl = GM, B 2 = 0. 


(B6) 
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One can then substitute this back into equations (Bl) 
and (B2) to hnd expressions for and 4/': 

GM 


d4> 

dr 

d^* 

dr 


(l + 5T)r2 
GM 


and 


(l + 5T)r2 
The PPN parameter is defined in 

7=1. 


(B7) 

(B 8 ) 


and one finds 


(B9) 


This shows that the definition of Newton’s constant (21) 
is indeed correct. 

Finally, we must check for the consistency of the as¬ 
sumption Bs = 0. Using equations (B4| and (B 6 ) we 
find that the solution for x is given by 


^2 = 


A2(r3A2 + 8GMTi/2) 

4r3(l-k5T) 


(BIO) 


at the stellar radius. This exactly matches the solution of 
equation (161 and therefore confirms that our assumption 
that B 3 = 0 is self-consistent. 


Appendix C: Behaviour of the Modified 
Lane-Emden equation near the Origin 


us begin by discussing the Lane-Emden equation. This 
is singular at the origin and so one typically starts the 
integration at some small number <5^1. This is fine for 
most ODEs, however the Lane-Emden equation possesses 
a special homology symmetry whereby several solutions 
are related through an appropriate transformation of ^ 
and 9 (see [S5M7] L We are not interested in the pre¬ 
cise form of these transformations here but one of their 
key implications is that the second boundary condition is 
redundant and the Lane-Emden equation may be recast 
as a first-order ODE. This means that by starting with 
initial conditions at some small 5 and not ^ = 0 one may 
evolve along a completely different class of solutions to 
the physical ones. In order to circumvent this problem, 
one can expand the Lane-Emden function near the origin 
as 


0(0 = i-ae+pe+o ie) 


(Cl) 


and one finds 


a = - and 
6 


(C2) 

(C3) 


The equation is then numerically solved by evaluating 
9{^ = d) and 9'{^ = 6) using equation (Cl) with these 
coefficients. 

Returning to the MLE, one finds 


In this Appendix we briefly discuss the numerical 
method of integrating the MLE. The MLE is a second- 
order ordinary differential equation (ODE) and so one 
may naively expect this to be a simple task whereby an 
explicit iterative numerical method such as the Runge- 
Kutta family of methods can be used with the two bound¬ 
ary conditions at ^ = 0 to construct the full solution. Let 


a = 

P = 


1 T 
6+4 
n(l -f I4T) 
120 


5nT2 

32 


(C4) 

(C5) 


In this work, we use these coefficients with d = 10 to 
solve the MLE numerically. 














